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The term "empirical predictor" refers to a two-stage predictor of 
a linear combination of fixed and random effects. In the first stage, 
a predictor is obtained but it involves unknown parameters; thus, 
in the second stage, the unknown parameters are replaced by their 
^H , estimators. In this paper, we consider mean squared errors (MSE) 

^^ ' of empirical predictors under a general setup, where ML or REML 

estimators are used for the second stage. We obtain second-order 
approximation to the MSE as well as an estimator of the MSE correct 
Cy ' to the same order. The general results are applied to mixed linear 

models to obtain a second-order approximation to the MSE of the 
empirical best linear unbiased predictor (EBLUP) of a linear mixed 
effect and an estimator of the MSE of EBLUP whose bias is correct 
,._^ , to second order. The general mixed linear model includes the mixed 

i]^ ' ANOVA model and the longitudinal model as special cases. 

in 

'^ \ 1. Introduction. We consider a general linear mixed model of the form 

O: (1.1) y = X(3 + Zv + e, 

'NT ■ 

^D I where y is an n x 1 vector of sample observations, X and Z are known 

c~| ■ matrices, /3 is a p x 1 vector of unknown parameters (fixed effects) and v 

j^ \ and e are distributed independently with means and covariance matrices 

^ ■ G and i?, respectively, depending on some unknown vector of parameters a. 

We assume that p is fixed and X is of full rank p {<n). Note that cov(y) = 

T. = R + ZGZ'. 

^ ' Problems involving multiple sources of random variation are often mod- 

H \ eled as special cases of (1.1). For example, in the well-known ANOVA 
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2 K. DAS, J. JIANG AND J. N. K. RAO 

model we partition Z as Z = (Zi, . . . , Zq) and v = {v[, . . . , v')' , where Zj 
is re X Tj, fj is Tj X 1, i = 1, . . . ,g, and vi, . . . ,Vq are mutually independent 
with Vi ^ N{0,ailri) and e ~ A^(0, o"o/n)- (For notational convenience we 
use (Tj rather than the customary af to denote the ith variance compo- 
nent.) Note that the AN OVA model is a special case of (1.1) with R = a^In, 
G = diag(a"ilri , . . . , CTqlr^) and a = (uo, o"i, . . . , cjg)'. The (dispersion) param- 
eter space under the ANOVA model is Q = {a:ai>0, i = 0,1,. . . ,q}. The 
well-known "longitudinal" random effects model [Laird and Ware (1982)] is 
also a special case of (1.1). In this case y = {y'l, . . . ,y[y with 

(1.2) yi = Xil3 + ZiVi + ei, i = l,...,t, 

where y^ is rej x 1, Xi is rii x p and Zj is rii x ri. It is assumed that the y^'s 
are independent, cov(t;j) = Gi, cov(ej) = Ri, where Gi and Ri depend on a, 
and Vi and Cj are independent. It follows that S = diag(Si, . . . ,St) with 
Sj = cov(yj) = T,i = Ri + ZiGiZ-. [Note that the longitudinal model (1.2) is 
not a special case of the ANOVA model and vice versa.] The well-known 
Fay-Herriot (1979) model widely used in small area estimation is a special 
case of the longitudinal model. The (dispersion) parameter space under the 
longitudinal model is Q = {a :T,i is nonnegative definite, i = 1, . . . ,t}. 

Estimation of linear combinations of /? and realized v from (1.1), say 
fi = I'P + m'v, for specified vectors of constants I and m is of considerable 
interest in many practical applications, for example, the estimation of quality 
index, longitudinal studies, the selection index in quantitative genetics, plant 
varietal trials and small area estimation [Robinson (1991)]. Henderson (1975) 
obtained the best linear unbiased predictor (BLUP) of fi under model (1.1) 
as 

t{a) =t{a,y) 

(1 3) 

= l'(3 + m'v = l'(3 + s{(T)'{y-XP), 

where 

/3 = /3(a) = {X'i:-^X)-^X'T.-^y 

is the generalized least squares estimator, or best linear unbiased estimator 
(BLUE), oi(5,v = v{cj) = GZ'J:-\y - X(5) and s{a) = Y.~^ZGm. 

The BLUP estimator (1.3) is unbiased in the sense of E[i((T, y) — /i] = 
under (1.1). The mean squared error (MSE) of t(o") is given by 

(1.4) MSE[t((7)] = E[t(a) - lif = gi{a) + g^ia), 

where 

g^{a) = m'{G - GZ'Y.-^ZG)m 
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and 

52(a) = [/ - X's{a)]'{X'T.~^Xr\l - X's{a)]- 

see Henderson (1975). Results (1.3) and (1.4) do not require normality of 
random effects v and e. 

The BLUP estimator t{a) depends on the dispersion parameters o", which 
are unknown in practice. It is therefore necessary to replace cr by a consistent 
estimator a to obtain a two-stage estimator or empirical BLUP (EBLUP) 
given by t((T). Methods of estimating a include maximum likelihood (ML) 
and restricted maximum likelihood (REML) under normality, the method 
of fitting-of-constants and minimum norm quadratic unbiased estimation 
(MINQUE) without the normality assumption; see Searle, Casella and Mc- 
Culloch (1992). The resulting estimators a are even and translation invari- 
ant, that is, d{y) = (j{—y) for all y and a(y + Xf3) = a{y) for all y and (3. 
Jiang (1996) proved that ML and REML estimators a obtained under nor- 
mality remain consistent without the normality assumption. 

Kackar and Harville (1981) showed that the EBLUP t{a) remains unbi- 
ased if a is even and translation invariant. This result holds provided that 
E[t((T)] is finite and v and e are symmetrically distributed (not necessarily 
normal). In particular, the two-stage estimator [5 = (3{a) is unbiased for (5. 

Kenward and Roger (1997) studied inference for the fixed effects, (5, un- 
der a general Gaussian linear mixed model y ~ N{Xf3, S) with a structured 
covariance matrix S = '^{a') depending on some parameter a. They used the 
REML estimator of /?, namely the two-stage estimator f3 = f3{a), where a 
is the REML estimator of a. A naive estimator of cov(/3) that ignores the 
variability in a is given by [X'T,~^ {a)X]~^ . Kenward and Roger (1997) de- 
rived a bias-adjusted estimator of cov(/3) and used it to derive a scaled Wald 
statistic, together with an F approximation to its distribution. The F ap- 
proximation performed well in simulations under a range of small sample 
settings. Kenward and Roger (1997) did not study the precise order of the 
bias of the adjusted estimator. 

Kackar and Harville (1981) studied approximation to the MSE of EBLUP 
t{a), assuming normality of the random effects v and errors e in the model (1.1) 
They showed that 

(1.5) MSE[t((T)] = MSE[i(a)] + E[t{a) - t{a)f 

for any even and translation invariant estimator a, provided that MSE[t(o")] 
is finite. It should be pointed out that, under very mild conditions, E[i((T)] 
and MSE[t((T)] are, in fact, finite [see Jiang (2000)]. It is customary among 
practitioners to ignore the variability associated with a and use the following 
naive estimator of MSE[i((T)]: 

(1.6) mse]<i[t{a)]= gi{a) + g2{a). 
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However, it follows from (1.4) and (1.5) that (1.6) can lead to significant un- 
derestimation. Therefore, it is practically important to obtain approximately 
unbiased estimators of MSE that reflect the true variability associated with 
the EBLUP estimators. This becomes particularly important when large 
sums of funds are involved. For example, over $7 billion dollars of funds 
are allocated annually on the basis of EBLUP estimators of school-age chil- 
dren in poverty at the county and school district levels [National Research 
Council (2000)]. 

Kackar and Harville (1984) gave an approximation to MSE[t((5")] under 
the general model (1.1), taking account of the variability in o", and proposed 
an estimator of MSE[i((T)] based on this approximation. However, the ap- 
proximation is somewhat heuristic, and the accuracy of the approximation 
and the associated estimator of MSE[i((T)] was not studied. Prasad and Rao 
(1990) studied the accuracy of a second-order approximation to MSE[i((T)] 
for two important special cases of the longitudinal model (1.2): (i) the well- 
known Fay-Herriot model (2.15) studied in Section 2.3 and (ii) the nested 
error linear regression model given by (1.2) with Zj = 1„. , a scalar Vi with 
vai{vi) = a\ and cov(ei) = aolm, i = I, ■ ■ ■ ,t. In the context of small area 
estimation n^ is the sample size in the ith area and t is the number of small 
areas. The nested error model may also be regarded as a special case of 
the AN OVA model with a single source of variation (g = 1), G = ailt and 
R = agin- Using the method of fitting-of-constants estimator a, Prasad and 
Rao (1990) showed that, for large t, 

(1.7) E[t{a)-t{a)]' = gs{a)+o{t~'), 

where 53(0") depends on cov((t). This leads to a second-order approximation 

(1.8) MSEa[t(<T)] =gi{a)+g2{a)+g3ia). 

The approximation is accurate to terms o{t~^), that is, the neglected terms 
are o(t~^). The 53(0") term is computationally simpler compared to an 
asymptotically equivalent term obtained from Kackar and Harville's approx- 
imation. Prasad and Rao (1990) also obtained an estimator of MSE[t((T)] 
given by 

(1.9) msepR[i(a)] =gi{a)+g2{a) + 2g3{a). 

The estimator (1.9) is approximately unbiased in the sense that its bias 
is o{t~^). Kackar and Harville (1984) proposed an alternative estimator given 
by 

(1.10) mseKK[t{a)] = giia) + g2ia) + g*3{a) 

for any even and translation invariant estimator a. The bias of (1.10) is 
0(i~^); that is, it is not approximately unbiased to terms o{t~^). 
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Harville and Jeske (1992) studied various MSE estimators under the 
AN OVA model with a single source of random variation (q = 1) and the 
REML estimator of a, including an estimator of the form (1.9). They re- 
ferred to the latter estimator as the Prasad-Rao (P-R) estimator. They 
appealed to Prasad-Rao's asymptotic results for its justification, but the 
latter results have been justified only for the special cases (i) and (ii). 
They also conducted a limited simulation study based on the simple one- 
way random effects model, yij = (3 + Vi + eij, j = 1, . . . , n^, i = 1,. . . ,t, us- 
ing a small balanced design (t = 9, rij = 2 for all i), a small unbalanced 
design (t = 9, ni = • • • = ng = 1, ng = 10) and a large unbalanced design 
(t = 21, ni = • • • = n2o = li "-21 = 50). The objective was to estimate the 
mean fi = P + vi. Simulation results indicated that the P-R estimator per- 
forms well when 7 = ai/ao is not small, but it can lead to substantial over- 
estimation for small values of 7 closer to the lower bound of 0, especially for 
small t (=9). Two partially Bayes estimators perform better than the P-R 
estimator when 7 is close to 0. 

Datta and Lahiri (2000) extended Prasad and Rao's (1990) results to the 
general longitudinal model (1.2) with covariance matrices Ri and Gi having 
linear structures of the form Ri = J2'j=Q(TjRijRij and Gi = Yl'j=Q'^jGijG[j, 
where uq = 1, Rij and Gij (i = 1, . . . , t; j = 0, 1, . . . , g) are known matrices 
with uniformly bounded elements such that Ri and Gi are positive definite 
matrices for i = l,...,t. They studied ML and REML estimators of a = 
((Ti, . . . ,(T|j)' and showed that an estimator of MSE[t((T)] of the form (1.9) is 
approximately unbiased when the REML estimator of a is used but not when 
the ML estimator is used. In the latter case an additional term involving the 
bias of the ML estimator a is needed for getting an approximately unbiased 
MSE estimator. Datta and Lahiri (2000) also gave explicit formulas under 
ML and REML for the two special cases (i) and (ii) studied by Prasad and 
Rao (1990). The underlying proof of Datta and Lahiri (2000), however, is 
not rigorous. 

The main purpose of our paper is to study the general linear mixed 
model (1.1) that covers the ANOVA model as well as the longitudinal model 
and derive a second-order approximation to MSE of EBLUP t{a) under 
REML and ML estimation of the variance components parameters a. We 
also derive approximately unbiased estimators of MSE[i((T)] and specify the 
precise order of the neglected terms. For ANOVA models with multiple 
sources of random variation, the components of a may have different con- 
vergence rates [Miller (1977) and Jiang (1996)]. As a result, rigorous proofs 
are quite technical and long. We have therefore only sketched the technical 
details in Section 5 of our paper, but the detailed proofs are available at the 
web site address given in Section 5. 

The remaining sections of the paper are organized as follows. In Sec- 
tion 2, we first present a general asymptotic representation of cj — a, where 
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a is obtained as a solution of "score" equations of the form dl{a)/da = 0, 
and a represents the true value of parameter vector. Normality assumption 
is not needed for this asymptotic representation. We then verify that the 
conditions underlying this representation are satisfied by solutions to the 
ML and REML score equations belonging to a parameter space O under the 
AN OVA model and normality. As another example, we show that the condi- 
tions are satisfied by the ML and REML estimators under the Fay-Herriot 
model and normality. In Section 3 we obtain a second-order approxima- 
tion to MSE[t(6")] under normality. The second-order approximation is then 
spelled out under the ANOVA model using ML and REML estimators a. 
We also verify the underlying key conditions for the special cases of the bal- 
anced ANOVA model and two special cases of the longitudinal model: the 
Fay-Herriot model and the nested error regression model. Section 4 gives an 
estimator of MSE[t((T)] correct to second order. The MSE estimator is then 
spelled out under the ANOVA model and the longitudinal model using ML 
and REML estimators a. Technical details are sketched in Section 5. 

2. Asymptotic representation of ct — cr. Throughout the rest of this 
paper, a represents the true parameter vector in places where there is no 
confusion; expressions such as dl{a)/da mean derivative with respect to 
a evaluated at a; in expressions such as E[5^^(a")/(9a"^], the expectation is 
taken at the true a, and the function inside E(-) is also evaluated at the 
true a. We first obtain an asymptotic representation of a — a, where a is 
obtained as a solution to "score" equations of the form 

(2.1) ^ = 

and then apply the general theory to the ANOVA model with REML and 
ML estimation of a. In Section 5.1 we sketch the proof of the asymptotic 
representation. Here l{a) may correspond to the logarithm of the restricted 
likelihood /r(o"), or the profile loglikelihood /p(cj) under model (1.1) and 
normality of v and e. 

Theorem 2.1. Suppose that: 

(i) l{a) =l{a,y) is three times continuously differentiable with respect 
to a = {ai,...,as)', where y = (yi, . . . ,y.„)'; 
(ii) the true a £ Q°, the interior of Q; 
(iii) 

(2.2) -oo<limsupAmax(^~^^^"^)<0, 

n— >oo 

where Amax means the largest eigenvalue, A = E[d'^l{a)/da'^] and D = diag((ii, 
. . . , dg) with di > 0, 1 <i < s, such that d* = mini<j<s dj — > oo as n ^ oo; 
and 
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(iv) the gth moments of the following are bounded (g > 0): 



1 

di 



dl{a) 



dai 



\/didj 



dH\ 



cr 



daj daj 



dn{a) 



d* 



k 

l<i,j,k<s, 

where Mijk{y) = sup^(zSg{a) \d^li^)/dai daj dak\ with Ssia) = {a : \ai - crj| < 
5d^/di, 1 < i < s} for some 5 > 0. Then (1) a a exists such that for any < 
p <1 there is a set B satisfying for large n and on B, a £ Q, dl{a)/d(T = 0, 
\D{a — a)\ < d^T'^ and 

(2.3) a = a-A~^a + r, 

where a = dl{a)/da and \r\ < d^ ^u^ with E(nf ) hounded, and (2) V{B'^) < 
cd^ ^^ , where r = (1/4) A (1 — p) and c is a constant. 

Note that Theorem 2.1 states that the solution to (2.1), o", exists and Ues in 
the parameter space with probabihty tending to 1 and gives the convergence 
rate of o" to o" as weh as the asymptotic representation (2.3), assuming that 
the true vector cr belongs to the interior of the parameter space G. 

2.1. REML estimation under the ANOVA model. The ANOVA model is 
given by 

(2.4) y = X(3 + Y,Z^Vi + e. 



The restricted loglikelihood under the ANOVA model with normality of v 
and e has the form 



(2.5) 
where cr 



Inicr) 



(l/2)[log(|r'Sr|)+y'Py], 



((To,(Ti, . . . ,o"g)', c is a constant, |T'ST| is the determinant of 
T'YiT, T is any n x {n — p) matrix such that rank(T) = n — p and T'X = 0, 



(2.6) 



P = T{T'T.TY^T' 



T.-~^X{X'T.-'^X)-'^X'T.-^, 



and S = Y11={] '^i^i with Vq = /„ and Vi = ZiZ^, i>l [e.g., Searle, Casella and 
McCulloch (1992), page 451]. The REML estimator of cr is a solution to (2.1) 
with I (a) = /r(o"). Section 5.2 sketches the proof that shows the conditions 
of Theorem 2.1 are satisfied, provided that the same conditions under which 
REML estimators are consistent are satisfied [Jiang (1996)]. The actual proof 
is somewhat lengthy and uses results on moments of quadratic forms in 
normal variables and di = \\Z[PZi\\2 for the ANOVA model, where ||-B||2 = 
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[tT{B'B)]^''^ for a matrix B. A quadratic form in y or n = y — Xf3 appears 
in the formulas for the first derivatives of /r, 

(2.7) ""' ^ 

= -[u'PViPu-tr{PV^)], 0<i<q. 

Note that n~ A^(0, S). Similarly, the second and third derivatives involve 
quadratic forms in n; see Section 5.2. 

2.2. ML estimation under the ANOVA model. The (unrestricted) log- 
likelihood has the form 

(2.8) lil3,a) = c-^[log{m) + iy-Xl3y^-\y-X(3)], 
where c is a constant. We have 

(2.9) ^^^ = r^-'y-X'j:-'Xl3, 

(2.10) ^^^ = hy- X/3)'S-V,S-i(y - XP) - tr(S- V,)], 

oai 2 

0<i<q. 

Solving dl{P, a)/df5 = for /3, we obtain from (2.9) /3(a) = {X'Y.-^X)~^X' x 
Y^~^y. Substituting /3(o") for j3 in (2.8), and using (2.6), we obtain the profile 
loglikelihood 

(2.11) h{a)=c-\[\ogm)+y'Py]. 

It now follows that the MLE of a is the solution of the equation 

Note that /p((t) is not a loglikelihood function, but Theorem 2.1 does not 
require l{a) to be a loglikelihood function, so we can take l{a) = Ip^cr) in 
Theorem 2.1. Section 5.3 sketches the proof that shows the conditions of 
Theorem 2.1 are satisfied with the same d, as in the REML case and the same 
set of conditions, provided p, the dimension of /3, is bounded as n increases. 
Again, quadratic forms appear in the formulas for the first derivatives of 
/p(a): 

^ = \[y'PV.Py-tr{^-'Vi)] 

(2.13) ""' 

= -[u'PViPu-ti(i:~^Vi)], 0<i<q. 

Similarly, the second and third derivatives involve quadratic forms in u] see 
Section 5.3. 
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2.3. The Fay-Herriot model. In Sections 2.1 and 2.2 we considered ML 
and REML estimations under the ANOVA model. Now we consider a differ- 
ent case, the Fay-Herriot model [Fay and Herriot (1979)]. This model has 
been considered by many authors; see Ghosh and Rao (1994) for a review 
and extensions. 

Suppose that yi is a scalar random variable such that 

(2.14) yi = XiP + Vi + ei, i = l,...,t, 

where the Uj's are i.i.d. ~ N(0,a), e^'s are independent such that Cj ~ 
N{0,(j)i) with known (pi, and fj's are independent of e^'s. Furthermore, Xi 
is a known p x 1 vector of covariates, and /3 is an unknown vector of re- 
gression coefficients. In the context of small area estimation, y^ denotes a 
survey estimate of the ith area mean /Xj and Cj denotes the sampling error 
with known sampling variance, var(ej) = (pi. Furthermore, /ij is modelled as 
fj-i = x[j3 + Vi with model errors Vi. 

Note that model (2.14) is not a special case of the ANOVA model. In fact, 
it may be considered as a special case of the longitudinal model introduced 
in Section 1. Model (2.14) may be written in matrix form as 

(2.15) y = Xf3 + v + e, 

where y = {yi, . . . , yt)', v = {vi,..., vt)' ~ A^(0, alt),e = (ei, . . . , e* )' ~ 7V(0, $) 
with $ = diag((/)i, . . . ^(pt), X \s t x p with ith row x^, and v, e are indepen- 
dent. 

Now consider REML and ML estimations under the Fay-Herriot model (2.15). 
In Section 5.6 we sketch the proofs that show the conditions of Theorem 2.1 
are satisfied if REML or ML estimators of a are used, provided that cr is 
positive and the (piS are bounded. 

3. MSE approximation. We now obtain a second-order approximation to 
the MSE of EBLUP t{a). Under normality the MSE of t{a) satisfies (1.5), 
that is, 

(3.1) MSE[t((T)] = MSE[t(a)] + ^[t{a) - t{a)f , 

where MSE[t((T)] is obtained from (1.4). It remains to approximate the last 
term on the right-hand side of (3.1). 

Theorem 3.1. Suppose that the conditions of Theorem 2.1 are satisfied. 
Furthermore, suppose that t{a) can be expressed as 

K 

(3.2) tia) = Y.Xkia)Wkiy), 

k=l 

where K = 0{d'^) for some a > 0, and the following terms are bounded for 
some 6 > 2 and 5 > 0: 
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(i) 

(ii) 
(iii) 

(iv) 



max E|H^fc(y)|' 

l<k<K 



max suplAipfcr) 



K 

E 

fc=l 



aAfc(cr) 



da 



K 

fc=ll^-'^l<5 



^^AfcC 



0- 



aa2 



where \\B\\ = [Amax(^'-B)]^/^ is t/ie spectral norm of a matrix B. If g > 
8(l + a)(l-2/6)"\ then 

(3.3) E[t{a) - t{a)f = E{h'A-^a)^ + o{d:^), 

where h = dt{a)/da, A = 'Ei[d'^l{a)/da], and a = dl{a) da . 

A sketch of the proof of Theorem 3.1 is given in Section 5.4. Note that 
normality is not required in Theorem 3.1. On the other hand, Theorem 3.1 
requires that the predictor t{a) have the form (3.2). In the next two sub- 
sections we show that this condition holds for balanced ANOVA models as 
well as for two longitudinal models. It is possible to replace (3.2) by some 
moment conditions on t{a) and its first and second derivatives, provided 
that one considers instead a truncated version of a, which is defined as a 
if |(t| < L„, and is a* otherwise with a* being a known vector in and L„ 
a positive number such that L„ — > 00 as n — > 00. The details of the latter 
approach are available at the web site given at the beginning of Section 5. 

3.1. ANOVA model. We first speU out E{h'A-^af for the ANOVA model 
with normality of v and e. We assume that the elements of the coefficient 
vectors / and m defining ^ = I' (3 + m'v are bounded, and \m\ = 0{\). In fact, 
m typically consists of only a finite number of elements equal to 1 and the 
rest equal to 0. For the balanced case /3 = {X'X)~^X'y does not depend 
on S. In this case h{a) =dt{a)/da = [V s{a)]'{I - Px)u = \y s{a)]'u+ lower 
order terms, where u = Zv + e, Vs{a) = ds{a)/da' and Px = X{X'X)~^X' . 
For the general unbalanced case /3 depends on S but the same expression still 
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holds, that is, dt{a)/da = \SJs{a)]'u+ lower order terms. Using Lemma 3.1 
below on higher moments of normal variables, we get 



¥.{h'A~^af = tr{[Vs(a)]'S[Vs((T)]A"^} + o(C^) 



(3-4) . . , . ._2 



when a is taken as dl-^{a)/da corresponding to REML. Note that d? is the 
diagonal element of the information matrix associated with ai and represents 
the "effective sample size" for estimating (Tj. 

Lemma 3.1. Let Ai and A2 be k x k symmetric m,atrices, and u ~ 
A^(0, S), where T, is k x k and positive definite. Then the following hold: 

(i) E[u{u'AjU - E{u'Aju)}u'] = 2^AjT., j = 1, 2. 
(ii) E[{u'Aiu - E{u'Aiu)}{u'A2U - E{u'A2u)}] = 2tr(AiSA2S). 
(iii) E[u{u'Aiu - E{u'Aiu)}{u'A2U - E{u'A2u)}u'] = 2tr(^iSyl2S)S + 
4EA1SA2S + 4SA2SA1S. 

(iv) Write Wj = X'-u, Wj = u'AjU, j = 1, . . . , s, where Xj and Aj are non- 
stochastic of order k x 1 and kxk, respectively. Then, for w = {wi, . . . , Wg)' 
and W = {Wi, . . . ,Wsy , 

E[w'{W-EW)f 

s s 

= tr(S^SvK) + 4^ ^ X'j1:{A,1:Ai + AiEAj)i:Xi, 

j=l 1=1 

where S^^ and S^y denote the covariance matrices of w and W , respectively. 

The proofs of (i)-(iv) are immediate from Lemmas A.l and A. 2 of Prasad 
and Rao (1990). 

It can be shown that (3.4) is also valid when a is taken as dl-p{a)/da 
corresponding to ML. Thus, using (1.4), (3.1), (3.3) and (3.4), we get 

(3.5) MSE[t(a)] = gi{a) + g2{a) + g^{a) + o(C'), 

valid for both the REML estimator and the ML estimator of a. 

Next, we show that the key condition (3.2) on t{cT) is satisfied for all bal- 
anced ANOVA models. Note that, based on the expression given by Propo- 
sition 3.1 below, all the other conditions of Theorem 3.1 are either trivial or 
known to be satisfied in the balanced case with normality [see Jiang (1996)]. 
Of course, the verification of (3.2) does not require normality. 

A balanced t(;-factor linear mixed ANOVA model may be expressed as 

(3.6) y = Xp + Y, Z,Vi + e, 

i&S 
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where X and Zj's have the fohowing expressions [e.g., Searle, Casella and 
McCuUoch (1992), Rao and Kleffe (1988)]: X = 0J1+/ 1^^'^ with (si, . . . , s^+i) e 
S^+i = {0, l}-+i, Zi = (8)^^+/ Ijl^ with (ii, . . . ,i^+i) G 5 C 5^+1, where ® 
denotes the operation of a Kronecker product, n; is the number of lev- 
els for factor Z, 1„ represents the n-dimensional vector of I's, 1^ = In and 
l\ = ln- The (w + l)st factor corresponds to "repetition within cells," and 
thus Syjj^i = 1 and iw+i = l,i £ S. We use for the element (0, . . . , 0) in 
Sw+i and let S = {0} U S. The covariance matrix of y then has the form 

S = aoln + ^ CTiZiZl 

ies 

(3.7) «.+! 

= E A.(g)j;',, 

where J^ represents the n x n matrix of I's, J^ = In and Jl^ = Jn, K = o^i if 
I E 5, and Aj = if i ^ 5. 

Searle and Henderson (1979) showed that S~^ has the same form, 

(3.8) E-i= ^ r,(g)J^^, 

where the coefficients Tj in (3.8) are computed by an algorithm. From a 
computational point of view, the Searle-Henderson algorithm is easy to 
operate. However, with such a form it may not be so easy to investigate 
theoretical properties of S~^, which are important to the current paper. 
Jiang (2004) gives an alternative derivation of (3.8) with explicit expres- 
sions for the Tj's; see Lemma 3.2. First, note that under the balanced model 
we have r^ = Ylii=Qni,i S 5. This allows us to extend the definition of r^ to 

all i E Suj+i- In particular, p = rs = ns;=o ^^ ^'^d n = tq = OS '^^ ^^ shall 
use the abbreviations {fc/ = 1} and so on for subsets of L = {1, . . . ,w + l}. For 
example, if A;, n G S^+i, then {ki = 1, ui = 0} means {I £ L:ki = l,ui = 0}. 
Also, for i, j £ Sw+i, j <i means ji <ii, 1 <l <w + 1. Finally, |-B| denotes 
the cardinality of set B. 

Lemma 3.2. For any balanced w-factor mixed ANOVA model (3.8) holds 
with 

(-,0, ^AJ V- (-i)lfa=i.i.=o}li(.<,) ] 

(3.9 Ti=[-]< 2^ — -— \, ieSy,+i. 

\nj [ .^^^ ao + J2kes mn/rk)l(^k<j) J 

Using Lemma 3.2, the following proposition can be proved. A sketch of 
the proof is given in Section 5.5. 
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Proposition 3.1. For any balanced mixed ANOVA model, the BLUP 
t{a) given by (1.3) can be expressed as (3.2) with K <1 + |S'|2"'+^ [hence 
K = 0{1) and the terms below (3.2) bounded]. 

It is known that (3.2) also holds for some unbalanced ANOVA models. 
For example, see the nested error regression model discussed in the next 
subsection, which is also a special case of the longitudinal model. 

3.2. The longitudinal model. For longitudinal models the spelled-out ex- 
pression for F,{h'A~^a)'^ in (3.3) [up to a term o{d~'^)] is given in Datta and 
Lahiri (2000). In the following we show that the key condition (3.2) in The- 
orem 3.1 is satisfied for two special (and important) classes of longitudinal 
models: the Fay-Herriot model and the nested error regression model. 

First consider the Fay-Herriot model (see Section 2.3). It is easy to show 
the following: 

t 
(3.10) l'P = Y,ai{a)yi, 

i=l 
t 







m'v = "y \miVi 






i=l 






t t 


(3.11) 




= ^h{(y)yi- X! bi,k{cr)yk, 

1=1 i,k=l 


where 








ai{a) -- 


ft / \ ^1 




V~t'^ + '^V ^^^i 




b^{a) -- 


- f ^ \ 




\cr + (j)iJ 




bi,k-- 


_„ ' '^ V' V ^^^'^ ' "" 




\a + <f>iJ \^-^<7 + 4>k) \cr + (t>k 



We have the following result. The proof is straightforward. 

Proposition 3.2. For the Fay-Herriot model (2.14), the BLUP t{a) 
given by (1.3) can be expressed as (3.2) with K = (t + l)'^ — 1 [and the terms 
below (3.2) bounded], provided that (i) the (pi's are bounded from above and 
away from 0; (ii) |xj|, I <i <t, are bounded, and so are [l[ and J2i=i l^-il; 
and (iii) limmi X^in{t~^Yll=i^i^i) > ^^ 

Next we consider the nested error regression model. Suppose that 
(3.12) yij =Po + x'ijf] + Vi + eij, j = 1, . . . ,ni] i = l, . . . ,t, 
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where j3 = {Pu)i<u<p~i and /3„, <u<p— 1, are unknown regression coef- 
ficients; Xij's are known vectors of covariates; Uj's are independent random 
effects witli distribution N{0,ai); Cjj's are independent errors with distri- 
bution A^(0,o"o), and v and e are independent. W.l.o.g. let iii > 1. For the 
nested error regression model (3.12) /'/? and m'v can, again, be expressed as 
(3.10) and (3.11), where 

6i(o-)=7nJ-i )1' , 



-1 



We have the following result. A sketch of the proof is given in Section 5.7. 

Proposition 3.3. For the nested error regression model (3.12) the BLUP 
t{a) given by (1.3) can be expressed as (3.2) with K = 0{t^) [and the terms 
below (3.2) bounded], provided that (i) (Tq > 0; (ii) p, ui, \xij\, 1 <i <t, 1 < 
j < ni, are bounded, and so are \l\ andJ2l=i |"ij|; (^nd (iii) liminf Amin(i~"'S'a) > 
0, a = 1, 2, where Si = E-=i ^i E]Li{xij-Xi)ixij-Xiy , S2 = E*=i E"=i ELi E^M' 
Xki){xij - XkiY and Xi = n^^ YJjLi Xij- 

Note that in both cases considered above d* can be chosen as \/t. 

4. Estimation of MSE. We now turn to the estimation of MSE[i((j)]. 
We obtain an estimator mse[i(6")] correct to second order in the sense of 
E{mse[t(a)]} = MSE[t((T)] + o{d-^). That is, the bias of mse[t{a)] in esti- 
mating MSE[t(o-)] is oid""^). 

First, we have from (3.1) and (3.3), 

MSE[i(a)] = MSE[t(a)] + E{h'A~^af + o(C^) 

= r]{a)+o{d^ ), say. 

We now define an estimator r} of 'r]{a) having the following property: 

(4.2) E(7?)=ry(a) + o(C'). 

It follows from (4.1) and (4.2) that the bias of fj in estimating MSE[t(o-)] 
is o(d~^). In addition to a and A defined in Section 2 (Theorem 2.1), 
let b = dr]{a)/da = (hi); B = d^rj{a)/da^ = (6,^); F = dH{a)/da'^, Hi = 
{d^l{a)/daida'^) and C = {a'A~^hi)i<i<s, where s is the dimension of a. 



MSE OF EMPIRICAL PREDICTOR 



15 



Also, let Q = D~^AD~^ and W 



\Wij') 



LetZ)~ia = (Ai),Z)-i/2(F- 



A)D ^/^ = (Ajj) and D ^HiD ^ = (Xijk)- We define the following vector, 
matrix and arrays: Uq = {ui), Ui = (uu), U2 = {ujm) and U3 = {ujkimn), 
where Ui = E(Aj), uu = E{XiXi), Ujti = E{XjkXi) and Ujkimn = E{XjkmXiXn). 
Note that ah of these are functions of a [e.g., A = A{a)]. The norm of an 
r-way array (r > 3) U, denoted by \\U\\, is defined as the maximum of the 
absolute values of its elements. (The norm of a matrix is defined in Theo- 
rem 3.1.) Define 

Ao{a) = -2b'A-^E{a), 

Ai{a)=b'A~^E{FA~^a), 

A2((t) = ^E{aA^^BA~^a), 



(4.3) 
(4.4) 
(4.5) 
(4.6) 
Finally, we define 



Asia) = -^b'A-^E{CA-^a). 



(4.7) f, = via)-J2A,ia), 

3=0 

provided that |f/| < cqc?^; otherwise, let fj = t?((T*), where cq and A are known 
positive constants and a* is a given point in 0. 



Theorem 4.1. The estimator fj given above satisfies the property (4.2) 
provided that 

(i) 7?(-) is three times continuously differentiable and the following are 
bounded: rj{a), \b\, \\B\\ and 

d'^rj{a) 



sup 

creSs(a) 



dai daj dak 



1 < ^5 j,k < s, 



where 6 is a positive number and Ss{(j) = {a : |o"j — (Tj| < 5d^/di, 1 <i < s}. 

(ii) The conditions of Theorem 2.1 hold with g > 8 + AX and l{a) four 
times continuously differentiable with respect to a. 

(iii) The gth moments of the following are bounded: 



1 



yjdjdk 



dH{a) 



and 



da.i daj dak 



E 



dH{a) 



, , , , sup 



dai daj da^ 
dH{d) 



1 



djdk 



d% 



cr 



dai daj dak 



daidajda^dai 



1 < i, j, A;, / < s. 



(iv) sup^g5_^(^) \\Q{a) - Q{a)\\ -^ and sup-^g^^^^^ \\Uj{a) - Uj{a) 
j = 1, 2, 3, as 6 ^ uniformly in n. 



0, 
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(v) |E(a)| is bounded and supg-go/^N |E(a)|(j=5- — E(a)| — > as 5 — > uni- 
formly in n. 

A sketch of the proof is given in Section 5.8. 

Corollary 4.1. // condition (v) of Theorem 4.1 is strengthened to 

(4.8) E(a)=0, 

then (4.2) holds, in which fj is given as in Theorem 4.1 except that in (4.7) 
the summation is from 1 to 3 . 

The estimator f) considered in Theorem 4.1 is truncated when its value 
exceeds some (large) threshold. Such a truncation is needed only for estab- 
lishing the asymptotic result. In practice one does not need to truncate the 
estimator (because it can be argued that for a given value of f/ there are al- 
ways constants cq and A such that \fj\ < cod^). On the other hand, a similar 
result may be obtained for rj without truncation, provided that a is replaced 
by its truncated version (defined above Section 3.1). The details of such a 
result are available at the web site given at the beginning of Section 5. 

4.1. REML and ML under the ANOVA model. In Section 5.9 we give 
sketches of a proof that shows the conditions of Theorem 4.1 are satisfied 
if the REML estimator of a is used, provided that the same conditions 
under which the REML estimators are consistent [Jiang (1996)] are satisfied. 
It can be shown that similar results also hold for ML estimation, but we 
omit the details. According to (3.5), in both REML and ML cases we have 
ri{a) = gi{a) +52(0') + gz{(y)- However, there is a difference between the two 
in terms of the spelled-out formulas for MSE estimation. This difference is 
due to the fact that (4.8) holds for REML but not for ML. 

First consider REML. Letting a = or, l{a) = /r(o") and A = A-^, we have 
E(aR) = 0, and the (i, j) element of Ar is given by —{l/2)ii{PViPVj). Fur- 
thermore, we have Ao(o") = 0, Ai((t) = 6'^r w^{a) + o{d~'^) , where w-^{a) = 

(w;o,R, • • -^Wq.K)' with Wi,B^ = -iT{A^[ix{PViPVjPVk)]Q<j±<q}, i = 0,...,q; 
^2(0") = —53(0") -|- o{d~'^), where 53(0") is given by (3.4) with A = A^ and 
Asia) = -b'A^^wnia) + oid-"^). Thus for REML we have Ej=o^i(^) = 
~93{(^) + o(d^^). It follows from (4.7) that for REML 77 = t/r, where 

(4.9) ?}r = 5'i('5"r) +52(<5-r) + 253 ((Tr), 

where itr is the REML estimator of a. The MSE estimator r/R given by (4.9) 
depends on the data y only through itr. An alternative MSE estimator that 
is data specific can be obtained by using the following estimator of 53(0"): 

(4.10) hi^R) = {y- Xpy[Vs{a)]'A^^[Vs{a)]{y - Xp)l=a^, 
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where (3 is the BLUE given below (1.3). It can be shown that E[^3((Tr)] = 
33(0") +o((i~^). The estimator (4.10) is obtained from (3.4) by replacing S 
hy t = {y-X(3^){y-X(3j,y. 

Now consider ML. For simplicity, we assume that rank(X) =p is bounded. 
Then, similarly, letting a = om, ^(ct) = /p(o") and A = Am for ML, we have 
E(aM) = -5'Mo(o-) = -fi'MO with gMQ,i = (1/2) tr[(i;"i - P)Vi], i = 0,...,q. 
Furthermore, Ao(o-) = 2b'A^gMo; Ai(o-) = b'A^WM - b'A^guo + o{d~'^), 
where wm is w^ with P replaced by S~^; A2(cr) = —gsia) + o{d~'^) and 
A3(a) = -b'A^^WM + o{d-^). Also, ^'^m^^mo = (5ffi/9c7)'A^^ffMo + o(C^) = 
510(0-) + o((i;:2), say. Thus for ML we have J2^j^o'^ji(^) = 5io(o-) - 93(0-) + 
o(d~^), hence 17 = f}M, where 

(4.11) ^M =5i('5-m) + 5'2('5-m) + 253(o-m) -5io(o-m), 

where o"m is the ML estimator of a. Similar to the REML case, a data-specific 
estimator can be obtained by using gsiau) instead of gsiau)- 

4.2. REML and ML under the longitudinal model. For longitudinal mod- 
els the spell-out of (4.7) [up to a term o(d~^)] is given by Datta and Lahiri 
(2000) for REML and ML estimation. Note that, similar to the previous sub- 
section, there is a difference between using REML and ML. In the following 
we give regularity conditions such that the conditions of Theorem 4.1 are 
satisfied for longitudinal models. The assumption that Gi and Ri are linear 
functions of a can be relaxed. 

We consider REML estimation. Similar results also hold for ML but we 
shall omit the details. Let a = (ctq, . . . ,o"g)'. Suppose that 

1. Gi and Ri are linear in a such that ||Gi||, ||-Ri||, \\dGi/daj\\, \\dRi/daj\\, 
< j < q, and ||S~ || are bounded, and, as a ^ a, maxi<j<i ||Gj(cr) — 
Gi{a)\\ -^ 0, \\Riia) - i?i(o-)|| -^ uniformly in t: 

2. o" G 0°, the interior of 0. 

3. The following are bounded: p, ni, \\Xi\\, \\Zi\\, \l\ and J2*i=i l"^i|- 

0, where Bi is the (g'-l- 1) x (q + l) matrix whose (j, k) element is tr(S^ Sj j x 
S~ Sj^fc) with Tiij = dTji/daj. Then the conditions of Theorem 4.1 are 
satisfied for the longitudinal model. A sketch of the proof is given in 
Section 5.10. 

5. Sketches of proofs. In this section we give very brief sketches of the 
proofs involved in this paper. These include proofs of the theorems and other 
technical details. The detailed proofs are available at the following web site 
address: http://anson.ucdavis.edu/~jiang/jp8.pdf. 
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5.1. Regarding Theorem 2.1. Let a* = o" — A~^a and B = BinB2, where 

Bi =m< {\x\/2)di-'} andi?2 = {M{i+\xnc\r+is'/V3)ai+\xnc\f < 

|A|/2} with e = D-^a, A = X^^^{D~^AD-^), r, = D-^{dH{a)/da'^ - A)D-^ 
and C = T^8iXij^k{Mijk{y)/didjdk} ■ It can be shown that, on B, l{a) — l{a^:) < 

if |L'(cr — cr*)| = 1. Since the function l{a) cannot attain its inaxiniuni 
over the set {a:\D{a — (T*)| < 1} at the boundary of the set, the maxi- 
mum must be attained in the interior. Thus, there is a solution to (2.1) 
in {a: \D{a — a^)\ < 1}. Let a be the solution to (2.1) closest to u*. It can be 
shown that, on B, a £ 0, dl{a)/da = and \D{a — a)\ < d*"''. Furthermore, 
by Taylor expansion it can be shown that, on ^B, a — a = —A~^a + ri + r2 
such that |ri| < d^ ~''ui, \r2\ < dZ ^U2 with E(n^), j = 1,2, bounded. Fi- 
nally, it can be shown that P(-Bi) < cid^ ~ and P(-B2) — C2d^^ for 
some constants ci and C2. 

5.2. Regarding Section 2.1. The following lemmas are used to verify that 
the conditions of Theorem 2.1 are satisfied. 

Lemma 5.1. Let Q be a symmetric matrix and ^ ^ N{0,I). Then for 
any g >2 there is a constant c that only depends on g such that 'E\S,'Q^ — 
EeQC\'^<c\\Q\\l 

Lemma 5.2. For any matrices A, B and C, we have \iT:{ABC)\ < \\B\\ ■ 
ll^lb • llC'lb; provided that the matrix product is well defined. 

Lemma 5.3. Let ai, hi he real numhers and 5i>Q such that ai<hi + Aq, 

1 <i < s, where A^ = Yfn=i^j0^i- If ^ = Yli]=i ^j < 1; then Oj < hi -|- (1 — 
Ay^Ab, l<i<s, where Ab = J2j=i^jbj ■ 

We also use the following expressions for second and third derivatives 
of /r((t): 

(5-1) ?^ = ltT{PV,PV,) - u'PV,PV,Pn, 

oaiOaj 2 

d'^Ri^) -u'PV,PV,PV,Pu 



dai daj dak 



(5.2) 



+ u'PVjPVkPViPu + u'PVkPViPVjPu 
- ^MPViPV,PVk)+tv{PV,PVkPVj)]. 



Note that the second and third derivatives are involved in the conditions of 
Theorem 2.1. 
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5.3. Regarding Section 2.2. In addition to (2.13), we have 

(5-3) P^ = ltv{J:~'V,J:-'Vj)-u'PV,PV,Pn, 

oaiOaj 2 

^^^^' u'PViPVjPVkPu 



(5.4) 



dat doj dak 

+ u'PVjPVkPViPu + u'PVkPViPVjPu 

-^[tr(S-V,S-V,S-V,.) 



+ tr(S-ViS-iT4S-V,)]. 
From these expressions we obtain the following relationships: 



dai daj dak dai daj da^ 2 

where h =tT{PViPVjPVk), h = tr{PViPVkPVj) and Jr is /^ with P re- 
placed by S~^, r = 1,2. We assume that p = rank(X) is bounded. Then it can 
beshownthat |tr(E-Vi)-tr(Pyi)| <pa^\\tr{J:-'^ViJ:~Wj)-tT{PViPVj)\ < 

2,p{a,aj)-^ and \tx{ll-^V^i:-%T.~^Vk) - tT{PViPVjPVk)\ <^p{^^<^i<yk)-^ ■ 
Thus, by the result of the previous subsection, it can be shown that the 
conditions of Theorem 2.1 are satisfied. 

5.4. Regarding Theorem 3.1. Let p = 3/4. By Theorem 2.1 and Taylor 
expansion it can be shown that t{a) — t{a) = —h'A~^a + r, where \r\ < 
d'^'^^u and ^{v?) is bounded. Thus, we have E[t(6-) - t{a)f = E(-)^1b + 
E(-)2lgc, where {-f denotes [t{a)-t{a)f. The first term = E(/i'A~^a)2lg + 

0(C^^^^''^) + 0{d~^P), while E(/i'^-ia)2lBc = 0{d-^-^) for some u>0. 

5.5. Regarding Proposition 3.1. First, the following identity can be es- 
tablished: 



(5.6) ^i -(P 



]xx'\i:-\ 
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By the definition of BLUP for v, Lemma 3.2 and (5.6), it can be shown that 
"^i = EfcGS„+i cnnWi^kV, where the r^'s are given by (3.9) and 

5.6. Regarding the Fay-Herriot model. For REML estimation, the re- 
stricted loghkehhood is given by Ir((t) = c— (1/2) (log |T'Sr| +y'Py), where 
c is a constant, T is as in Section 2.1 and P = T{T'T,T)~^T' = the middle 
term of (5.6) with S = 0"/ + $. Suppose that o" > and the (j)iS are bounded. 
Then it can be shown that the conditions of Theorem 2.1 are satisfied with 
D = d = \ft. A similar result can be proved for ML estimation, in which case 
one considers the profile loghkehhood Zp((t) = c — (1/2) [log |S| + y'Py]. 

5.7. Regarding Proposition 3.3. First note that X; = (l„.Xi), where the 
jth row of Xi is x-^. Also, we have Sj = cro/n, + o"i Jn^, thus S^"^ = Aj~^/„. + 
7njA,"^(/„, -J„J, where 7 = 0-1 /ctq, \i = \i{a) = aQ + niai and Jn, = JnJni. 
Therefore, we can write 

Ev'v-i V — ( 
.^^^^^* ^'-\B C + ^D 

where A = Yd^-^ Ui/Xi, B = J2J^^ x-ln^/Aj, C = J2i=i x'^Xi/Xi and D = Y!i=i{ni/Xi)x[{In^ ■ 
Jn,)xi. Thus, 

(5.7) (i:^;5:-^,)"=(_Ab "™' 

w\ieie Q= [A -B'{C + -iD)-^B]-^ and R = {C + -f D - A'^ B B')-\ It can 
be shown that AC — BB' > 52/2A|^, where Am = maxjAj = (Tq + nmaxO-i 
with rimax = maxj n j . It follows, by conditions (ii) and (iii), that ||i?|| < 
AM/t(5i7 + 62), where 5a, a = 1,2, are some positive constants. Then, using 
the identity Q = A"^ + A~^B'RB, one can show \\Q\\ < c^Xu/t), where c is 
a constant. 

5.8. Regarding Theorem 4.1. The proof of Theorem 4.1 requires the fol- 
lowing lemmas [see Jiang, Lahiri and Wan (2002)]. 

Lemma 5.4. For any nonsingular matrices P, Q and nonnegative inte- 
ger q, 
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Lemma 5.5. Let P, Q be matrices such that P is nonsingular and \\Q ■ 
P\\ < (3||P-i||)-i. Then Q is nonsingular and \\Q~'^\\ < {3/V2)\\P~^\\. 

Let A be the set such that the fonowing hold: 



1 



dH\ 



o) 



Vdjdk 



\/didj 
dn{a 



dui doj 



E 



dH{a) 



dui doj 



dai daj dak 



E 



dH{ 



(T 



dai daj dau 



<dl, l<i,j<s, 

<dl, l<i,j,k<s. 



Let £ = Ar\B. Let p = 3/4 in Theorem 2.1. It can be shown that PiS") < 
cd^^^ , where r = 1/4. By Taylor expansion, it can be shown that the fol- 
lowing holds on £: 

r]{a) = r]{a) - 2b'A-^a + b'A'^fA~^a 

+ ^[a'A-^BA-^a - b'A-^CA-^a] + r, 

where \r\ < d* ^u and E(ti) is bounded. Thus it can be shown that ^r]{a)\e = 
r]{a) + J2j=o^j{<^) + o{d~^). On the other hand, we have the following ex- 
pressions: 

1 



Ao(o-) = -2^—bi{a)wij{a)uj{a) 



*j 



di 



1 



"'<^'^.£*^ 



bi{a)wij{a)wki{a)ujkiia), 



A2(o-) = - Y^ -T-rbjk{a)wij{a)wki{a)uii{a), 
^ i,j,k,i ^i"'^ 

1 1 

^3(0-) = -7^ X! -T-rbi{<y)Wij{a)wkl{a)w,nn{(r)Ujklmn{(r)- 
i,j,k,l,m,n -^ 

With these it can be shown that EAj{a)l£ = Aj(cj) + o{d~'^). It follows that 
Ef/ls = i]{a) + o{d^'^). Finally, we have E|f/|l£:c = o(d~^). 



5.9. Regarding Section 4.1. It suffices to show that (3) and (4) hold. Let 
i' , j' , k' be a permutation of i, j, k and w.l.o.g. let di' = di' A dj' A dk' = 
di A dj A dk- Then by Lemmas 5.1 and 5.2, it can be shown that 



E 



1 



^Jdjdk 



[y'PVePV.'PVk'Py - ^{y'PV,,PVj,PVk'Py)] 



<c. 
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Here c represents a constant whose value may be different at different places. 
Similarly, it can be shown that | tr(PViiPVjiPVk')\ < c{didjdk/di V dj V d^), 



didjdkdi ^(zsg 



e( -nrr-T ^up \y' PVi> PVj> PVk> PVi> Py\ ] < c 



and {dl/didjdkdi)\tiiPVi,PVj,PVk>PVi,)\<c. 

As for (4), first note that P{a) = -{l/2){tr{HGiHGj)/didj)o<ij<g. It 
can be shown that 

I ti{HGiHGj) - tv{HGiHGj)\ 

< E \^k - ak\MHGkHG^HGj)\ 

k 

+ Y.\^i-''i\\^'i^Gi^GjHGi)\ 

I 

+ J2\ak-ak\\al-al\\tv{HGkHG^HGlHGj)\, 

k,i 

\tr{HGkHGiHGj)\ <2a^^d^dj and \tv{HGkHGiHGiHGj)\ < 4(0-^0-;)"^ x 
didj. Thus supa-g^^ ||P(cr) - P{a)\\ ^ as 5 ^ 0. Note that Ui{a) = -P{a). 
Similarly, one can show that sup^.^^^ ||f^2(5') — [/2(<7)|| ^ as 5 ^ 0. Finally, 
it can be shown that 

""iiklm = ^r ,\ , I MPV^PVjPVl)+tr{PViPVlPVj)]tT{PVkPVm) 

Zdjakaiam I 



+ E tliPVaPVbPVcPVkPVm) 
a,b,c 

+ E triPVaPVbPV,PVmPVk)\, 
a,b,c J 



where the summation is over all a, b, c which is a permutation of i, j, I. It 
follows that sup5.g5^ llf^slc^) — ^3(o")|| — > as 5 — > 0. 

5.10. Regarding Section 4.2. First note that the formulas derived in Sec- 
tions 2.1 and 5.2 for /R(cr) and its derivatives hold for the general linear 
mixed model (1.1), which includes the longitudinal model. Next, note that 
the matrix P of (2.6) can be expressed as P = S~^ + A, where ||A||2 is 
bounded. These results and Lemma 5.1 are used to verify the moment con- 
ditions involved. 
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As for the conditions regarding r/ and its derivatives, we have the following 
expressions: r]{a) = gi{a) + 52(c) + (73(0"), where 



i=\ 






, i=l 



i=\ 



93(fT) = ^tr| 



„ (S- ZjGjmj 
da ' 



d_ 
da 



'llT^ZiGimi) 



With these expressions one can verify the conditions regarding 77 b and ||-B||. 
Finally, for condition (iv) of Theorem 4.1 we have, for example, Q — Q = 
-(l/2t)[tr(PSjPSfc)]o<j,fc<g, where Sj = dT./daj. Note that P-P = P{E- 
T,)P, where P is P with a replaced by a and so on. 
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